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The need to smoothly cover a computational domain of interest generically requires the adoption of 
several grids. To solve the problem of interest under this grid-structure one must ensure the suitable 
transfer of information among the different grids involved. In this work we discuss a technique that 
allows one to construct finite difference schemes of arbitrary high order which are guaranteed to 
satisfy linear numerical and strict stability. The technique relies on the use of difference operators 
satisfying summation by parts and penalty techniques to transfer information between the grids. 
lO ' This allows the derivation of semidiscrete energy estimates for problems admitting such estimates 

, at the continuum. We analyze several aspects of this technique when used in conjuction with high 

. order schemes and illustrate its use in one, two and three dimensional numerical relativity model 

problems with non-trivial topologies, including truly spherical black hole excision. 

I. INTRODUCTION 

Many systems of interest have a non-trivial natural topology that a single cubical computational domain cannot 
^ ■ accommodate in a smooth manner. Examples of these topologies in three-dimensional settings are x R (or a subset 
' of it), encountered when dealing with spacetimes with smooth outer boundary and inner boundaries where required 
I to excise singularities, and 5^ x or S*"^ topologies, commonly found in cosmological problems. 
' The need to treat these scenarios naturally leads one to consider multiple coordinate patches in order to cover 
the region of integration. These, in turn, translate into having to adopt multiple grids at the implementation level. 
Each one of these represents a region of discrete space, a patch, which might come equipped with a discrete Cartesian 
coordinate system or discrete charts; that is, invertible maps from discrete space to regions of Z'^ (i.e. the integers 
labeling coordinate grid points in each direction). 
^ Covering a spacetime by charts is commonly done at the continuum level when considering the differential geometry 

I of the spacetime (see for instance Q). These charts are usually thought of as defining a map of a portion of the 
^jj^ spacetime into a subset of , and the combination of charts (which usually overlap in some regions) covers the whole 
• . , spacetime. Points belonging to an overlapping region are considered as belonging to any of the involved charts. Here 
a well defined coordinate transformation between the charts is naturally defined by the combination of the maps and 
their inverse in between the spacetime and the charts. 
; ^ • At the discrete level one can in principle adopt an analog of the above construction. However, it is often the 
Ci ' case that in the overlapping region a grid point in one of the charts does not have a corresponding one in the 
other. Consequently, the coordinate transformation is not defined. This presents a problem at the practical level 
as communication between patches must take place. This issue is commonly solved in two different ways: (I) By 
introducing further points via interpolation where needed, (H) By considering patches that only abut, i.e. do not 
overlap. 

In the first case - commonly referred to as overlapping-grids approach - non-existent points in a given grid within the 
overlapping region are defined where needed by appropriate interpolations. Although this can be done in a straight- 
forward manner with a relatively simple multiple-grid structure, the drawback of this approach is the introduction 
of a new ingredient - the interpolation - which does not have a counterpart at the continuum. This complicates 
the assessment of stability of even simple evolution problems as the details of the interpolation itself arc intertwined 
with any attempt in this direction in an involved manner. As a consequence, there exist few stability proofs for such 
evolution schemes and have so far been restricted to one-dimensional settings. Notwithstanding this point a number of 
implementations in Numerical Relativity, where the possible truncation-error driven inconsistencies at the interfaces 
are dealt with by introducing a certain amount of dissipation or filtering, make use of this approach with good results 
(see for instance, HHIIQ). 

In the second case - commonly referred to as multi-block approach - grids are defined in a way such that there is 
no overlap and only grid points at boundaries are common to different grids. This requirement translates into having 
to define the multiple grids with greater care than in the previous option. This extra effort, however, has as one 
pay off that schemes preserving important continuum properties can be constructed. In particular, this allows the 
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construction of stability analyses which are similar to those of a single grid. More explicitly, following Abarbanel, 
Carpenter, Nordstrom, and Gottheb 0,0,13 one can construct schemes of arbitrary high order for which semi-discrete 
energy estimates are straightforward to derive in a general way. The availability of stability results for this second 
approach makes it a very attractive option in involved problems -like those typically found when evolving Einstein 
equations- where schemes eliminating spurious sources of instabilities provide a strong starting point for a stable 
implementation of the problem. 

In this paper we discuss and analyze the use of this multi-block approach in the context of Numerical Relativity. At 
the core of the technique to treat outer and patch interfaces is the addition of suitable penalty terms to the evolution 
equations 0, @- t;he case of hyperbolic systems these terms penalize the possible mismatches between the 
different values the characteristic fields take at the interface between several patches. 

Not only does this method provide a consistent way to communicate information between the different patches but, 
more importantly, does so in a way which allows for the derivation of energy estimates at the semi-discrete level. 
Consequently, numerical stability can be ensured for a large set of problems. These estimates can be obtained with 
difference operators of any accuracy order, provided they satisfy the summation by parts (SBP) property and the 
penalty terms are constructed appropriately. 

In this work we discuss this technique in a context relevant to numerical relativity, analyze its properties and 
illustrate it in specific examples. In particular we show results for the case of the x R topology used in black hole 
excision techniques. 

This work is organized as follows. Section II includes a description of the numerical analysis needed to attain 
stability in the presence of multiple grids and summarize how the penalty method of Refs. 0, 0, allows for 
achieving this goal. 

In Section III we study some aspects of Strand's high order operators satisfying SBP with respect to diagonal 
norms, when combined with the penalty technique. We find that in some cases, typically used operators that minimize 
the bandwidth have a very large spectral radius, with corresponding limitations in the Courant-Friedrich-Levy (CFL) 
factor when used in evolution equations. We therefore construct operators that minimize the spectral radius instead. 
Additionally, we examine the behavior of the convergence rate and the propagation behavior that different modes 
have when employing different higher order operators. 

In Section IV we present and analyze different tests relevant to numerical relativity employing derivative operators 
of different order of accuracy and the penalty technique to deal with multiple grids. These tests cover from linearized 
Einstein equations (in effectively one-dimensional scenarios) to propagation of three-dimensional fields in curved 
backgrounds. 

We defer to appendices the discussion of several issues. Appendix A presents details of the higher derivative 
operators and diagonal norms which we employ in this work. Appendix B discusses our construction of high order 
dissipative operators which are negative definite with respect to the corresponding SBP scalar product. Last, Appendix 
C lists some useful properties that finite difference derivative operator satisfy, which help in our construction of 
dissipative operators. 

II. INTERFACE TREATMENT FOR SYMMETRIC HYPERBOLIC PROBLEMS IN MULTIPLE 

BLOCKS 

As mentioned, we are interested in setting up a computational domain which consists of several grids which just abut. 
This domain provides the basic arena on which symmetric hyperbolic systems are to be numerically implemented. 
The basic strategy is to discretize the equations at each individual grid or block, treating boundary points in a suitable 
way. Boundary points at each grid either represent true boundary ones from a global perspective or lie at the interface 
between grids. In the latter case, since these points are common to more than one grid the solution at them can be 
regarded as multi-valued. As we show below, this issue can be dealt with consistently and stably, ensuring that any 
possible mismatch converges to zero with resolution. 

At the core of the technique is the appropriate communication of these possibly different values of the solution 
at the interfaces. Intuitively, since we are dealing with symmetric hyperbolic systems, a natural approach would 
be to communicate the characteristic variables from one domain to the other one. However, this is not known to 
be numerically stable. There exists nonetheless a technique based on this strategy which does guarantee numerical 
stability 0,013. This relies in adding penalty terms to the evolution equations of characteristic fields which penalize: 
a) in the interface case the mismatch between the different values each characteristic field takes at the interface of 
several grids; b) in the outer boundary case the difference between each incoming characteristic field and the boundary 
conditions one wishes to impose to it. 

These penalty terms are constructed so as to guarantee the stability of the whole composite grid if it can be 
guaranteed at each individual grid through the energy method. To this end, the use of schemes with difference 
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operators satisfying SBP are employed. Hence, on each single grid there exists a family of natural semidiscrete 
energies, defined by both a symmetrizer of the continuum equations and a discrete scalar product with respect to 
which SBP holds \M- One can then define an energy for the whole domain by simply adding the different energies of 
each grid. The use of operators satisfying SBP allows one to get an energy estimate, up to outer boundary and interface 
terms left after SBP. The penalties allow to control their contribution, thus obtaining an estimate for the global grid. 
This is achieved if the contribution to the time derivative of the energy due to the interface and outer boundary terms 
(in the latter case when, say, homogenous maximally dissipative outer boundary conditions are imposed) left after 
SBP is non-positive. When these terms are exactly zero the penalty treatment of Refs. ^0,13 is, in a precise sense, 
"energy non-dissipative" . On the other hand, if these terms are negative the scheme is numerically stable but at fixed 
resolution a damping of the energy (with respect to the growth one would obtain in the absence of an interface) in 
time arises. This damping is proportional to either the mismatch of a given characteristic variable at each interface 
or its failure to satisfy an outer boundary condition. As we describe below, these interface and boundary terms left 
after SBP are controlled precisely by the mentioned penalties, each of which depends on: the possible mismatch; the 
characteristic speeds, the corresponding SBP scalar product at the interface, the resolution at each intervening grid, 
and a free parameter which regulates the strength of the penalties. 

Next, we explicitly describe how this penalty technique allows one to derive semidiscrete energy estimates. We 
first discuss in detail the one-dimensional example of an advection equation on a domain with a single interface. The 
more general case of systems of equations in several dimensions follows essentially the same principles, applying the 
Id treatment to each characteristic field. We illustrate this by discussing a general constant-coefficient system in a 
given two-dimensional setting. From this, the generalization to the three-dimensional general case is straightforward, 
and we therefore only highlight its salient features. 



A. A one-dimensional example 

Consider a computational domain represented by a discrete grid consisting of points i = imin ■ ■ ■ imax and gridspacing 
h covering a: £ [a, 5]. A Id difference operator D on such a domain is said to satisfy SBP with respect to a positive 
definite scalar product (defined by its coefficients ct^ ) 

{u,v) = h'^UiVjUij , (1) 

if the property 

{u,Dv) + {v,Du) = {uv) \l 

holds for all gridfunctions w, v. The scalar product/norm is said to be diagonal if Uij = crM(5ij j37|. One advantage of 
Id difference operators satisfying SBP wrt diagonal norms is that SBP is guaranteed to hold in several dimensions 
if the Id operator is used on each direction (which is not known to hold in the non-diagonal case in general) |lClj| . 
Even in Id, in the variable coefficients and non-diagonal case the commutator between D and the principal part 
might not be bounded for all resolutions (something that is generically needed for an energy estimate to hold) 
pT]. Another advantage is that the operators are, for a given order in the interior, simpler in their expressions. The 
disadvantage is that their order at and close to boundaries is half that one in the interior, while in the non-diagonal 
case the operators loose only one order with respect to the interior |^ IT^ . Throughout this paper we will mostly 
restrict our treatment to the use of diagonal norms. 

As an example of how to impose interface or outer boundary conditions through penalty terms, we concentrate 
next on the advection equation for u propagating with speed A, 

u = AdxU . (2) 



1. A domain with an interface 

Consider the interval (—00,00) with appropriate fall-off conditions at infinity. We consider two grids: a left one 
covering (— oo,0], and a right one covering [0, 00). We refer to the gridfunction u on each grid by and (cor- 
responding to the left and right grids, respectively). Both of these gridfunctions have a point defined at the a; = 
interface and they need not coincide there, except at the initial time. Therefore, the numerical solution will in principle 
be multivalued at a; = 0, though, as we shall see, the penalty technique is designed to keep this difference small. 
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The problem is discretized using on the right and left grids, respectively, with gridspacings h^h^ -not necessarily 
equal- and difference operators Z?' , satisfying SBP with respect to scalar products given by the weights cr' , at 
their individual grids. That is, these scalar products are defined through 
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The semidiscrete equations are written as 



u\ = KD'u\+^-^{ul-u',), (3) 

u\ = Ai?X + ^K-^S)- (4) 
n ctqq 

Notice in the above equations the second term on each right hand side, which constitutes the penalty added to the 
problem. They are defined by the possible mismatch, the grid-spacing, the inner product employed and the free 
parameters {S\ S"^}, which will be determined by requiring an energy estimate to hold. 

We define a natural energy for the whole domain, which is the sum of the energies for each grid (for this simple 
example with a trivial symmetrizer). 

Taking a time derivative of this energy, using the semidiscrete equations (|3I4|I and the SBP property one gets 

Et = {K- 2S'){u'^f + (-A - 25'')K)2 + 2(5' + S')u'^ul . (5) 

In order to get an energy estimate the above interface term (i.e., the right-hand side of Eq. |SJl) must be non-positive 
for all Wq,Uq. It is straightforward to check that this is equivalent to the three following conditions holding: 

A -25; < (6) 
-A-25,r < (7) 
(A + 5^-5;)^ < (8) 

From there, it is clear that we need A + 5r — 5; =0. And with this condition the other two become 5; + 5^ > 0. 
There are three possibilities: 



• Positive A.- we can take 



Si = K + 5, 5, = (5, with ^ > (9) 

The time derivative of the energy with this choice becomes 

Et^-{v}o-ulf{K + 25)<Q 

• Negative K: similarly, we can take 

5^ = -A + (5, Si = 5, with (5 > ^ (10) 

The corresponding time derivative of the energy with this choice becomes 

Et^{u''a--ulf{K~25)<{) 

• Vanishing A: this can be seen as the limiting case of any of the above two, and we can take 

Si = Sr = (5, with 5>Q . (11) 

Hence, 

Et^^{u'^-ulf26<{) 
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The coefBcients Si and Sr need not be equal, but the following symmetry is important, under the change A — > —A 
we should have Si — *■ 5*^ and vice-versa, since it transforms incoming modes to outgoing ones. This is clearly satisfied 
by the choice above. 

Smumarizing, there is a freedom in the penalty factors of Eas. l|3l4l) . encoded in the parameter S, which has to satisfy 
6 > — |A|/2. If (5 = — |A/2| there is no interface term in the estimate (that is, exact energy conservation in the above 
model), while if If J > — |A/2| there is a negative definite interface term in the estimate (which represents a damping 
in the energy proportional to the mismatch). 

As we will see below, one proceeds similarly in the more general case of systems of equations in several dimensions. 
The penalty terms are applied to the evolution equation of each characteristic mode, with factors given by Eas. H9llUlllf) 
(where A in the general case is the corresponding characteristic speed). 

2. A domain with an outer boundary 

The penalty method also allows to treat outer boundaries in a similar way. As an example, consider again the 
advective equation Eq.(|21), but now on the domain (— cxd, 0], with gridpoints i = — oo ... 0. Assume A > 0; boundary 
conditions therefore need to be given at x = 0; say u{x = 0,t) = g. The semidiscrete equations are written as 

S T 

ii, = ADu, + ^^{g-uo). (12) 
tiaoo 

(13) 

Defining the energy to be 

E = 

its time derivative is 

E = {A-2T)ul + 2guoT, (14) 
< iA-T)ul + Tg\ (15) 

As in the interface case with positive speed [c.f. Eq.(|5J], we can therefore take T = A + 6. For the homogeneous case, 
(7 = 0, the equality (|14l) holds and we have 

= (A - 2T)ul , 

indicating that for 6 > —A/2 the energy will not increase. For the non-homogeneous case {g ^ 0) the inequality IjlSf) 
yields 

E < A.g2 + 5{g^ - u^) ; 

Note that for 5 = one trivially recovers the continuum estimate. For other values of S the consistency with the 
continuum estimate follows from the observation that uq converges to g to the s + 1-th order if the SBP derivative 
operator has accuracy s at the boundary point J^]. This implies that both the numerical implementation and the 
corresponding energy estimate are consistent with those defined at the continuum level. 

At this point we find it important to remark the following. Notice that just having an energy inequality is not 
enough, as one further needs to ensure consistency of the discrete equations with respect to the continuum ones. In 
the case of the penalty approach this not straightforward as the penalty term diverges when the grid size decreases 
unless u converges to g sufficiently fast (as mentioned above, u does converge to g fast enough if u is an incoming 
mode). In fact, the penalty term can be viewed at the continuum as approximating the original equation and boundary 
conditions through the introduction of a suitable delta function at the boundary. The argument of the delta function 
must be consistent with the underlying problem. For instance, if one were to put a penalty term on a boundary where 
the mode is outgoing -and thus the value of the function there is determined by the evolution itself- the inconsistency 
would manifest itself through a lack of convergence. In the above example this would be the case if we take A < 
but insist in putting a penalty term at a: = 0. Note that in such situation the energy inequality would still hold if 
T > A/2. This would imply that the numerical solution is still bounded in the norm, but no more than that; 
indeed, numerical experiments show that a high frequency solution traveling in the incoming direction (that is, with 
velocity opposite to that one at the continuum, see Section [ill C|) . whose amplitude depends on the size of the penalty 
term, is generated. Naturally, if T < A/2 the energy inequality is violated and the solution blows up exponentially 
with a rate increasing with the highest frequency that can be accommodated by the grid being employed. 



6 



B. The two-dimensional case 



Consider now the system of equations 

ii = A^'^^,u = A'd^u + AydyU , 

where u is a vector- valued function, A^^A^ are symmetric and, for simpUcity, constant coefficient matrices. The 
domain is composed of two grids: a left one covering (—00, 0] x (—00, 0] and a right one covering [0, 00) x {—00, 0], with 
an interface at x = and an outer boundary dX y — (see Figure At y = we impose homogeneous maximally 
dissipative boundary conditions by setting to zero the incoming characteristic fields. 
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FIG. 1: Example of a multi-block domain in two dimensions. Only the spacing in the vertical direction need be the same in 
both grids so as to ensure boundary points -represented by grey hexagons- coincide. 



We assume that the scalar products associated with the different Id difference operators are diagonal. This, as 
already mentioned, ensures certain properties that guarantee an energy estimate. 

We denote the gridspacings, difference operators and associated scalar products corresponding to the x, y directions 
by hx,hy, Dx, Dy, and cr^^), '^{y)y respectively. These quantities need not coincide on the different subdomains, except 
for the gridspacing in the transversal direction at an interface (so that gridpoints belonging to different grids align 
with each other). To all quantities we add an Z or r supraindex to denote quantities belonging to the left or right 
domain. Thus, the only condition we require is hy — hy =: hy. 

On each subdomain the 2d scalar product is defined as the product of the scalar product on each direction, 

{u, v) = hxhy'^{u^j,Vi■j)(T^x)i'7^y)j (16) 

where [u, v) is the pointwise Euclidean scalar product of two vectors. 

As in Id, the total energy is defined as the sum of the energies on each subdomain (in this case with a trivial 
symmetrizer, as the equations are already in symmetric form): 



where 



E = E' 



i<0 j<0 

E"' = Khy^^iKj,ulj)al^^^<7ly^^ . (18) 

i>0 j<0 



The evolution equations are a combination of Eas. H3l4|l and Ea. (|I2|l . 



= A'^Dlui^ + T|^^(-S. - O - Tr^^<o (19) 

- -0,) - J^T^^o (20) 
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In the above expressions, 5*', S^jT^T"^ are operators (as opposed to scalars), since we are dealing with a system of 
equations. The first two correspond to penalty terms added to handle grid interfaces while the latter two for imposing 
outer boundary conditions. The goal of these operators is to transform to characteristic variables and apply to the 
evolution equation of each characteristic mode suitable penalty terms, as in Eas. (|3l4ll2|l . 

Taking a time derivative of the energies defined in Ea. (|17ll8|l . using the evolution equations H19I20|1 . and employing 
the SBP property along each direction, one gets 

= hyY, a[y^^ [(4^^, [A^ - 2S'W,^^) + 2iu',^^, S%^)] + h^Y^ <jUM,o^ i^' " 2T'H,o) (21) 

j<0 i<0 

- /^yE^feb- [K,,(-^"-25'-)uS,,) + 2K,,5'-<,)] +/i;E^W^Ko:(^"-2^^^^ (22) 

j"<0 i>0 

where we have assumed that S and T are hermitian matrices. 
In order to control the interface terms in + E"^ we can take 

(J, ^- 

{y)j 



= ^ [(-A- + S-)P^ + 5+ PI + <50Po] ; 
{y)j 

where a sum over the index a is assumed, and {P^,P°,Po} are projectors to the sub-spaces of eigenvectors of 
with eigenvalues {A+, A^^, A"} respectively. With this choice E^ + E"^ becomes 



E^+E'' = h, 



'E[(A.- 



2<5; 



a,l a,r \ \ 2 



lP-(A^ + 25+ 



]<0 



'■i.Ol 



(23) 
(24) 



i<0 i>0 

Clearly, in order to obtain an estimate, the following conditions must be satisfied, 

A--2^-<0, A+ + 2(5+>0, (5°>0; 

which is analogous to the one-dimensional case, Eas. H9ll0lll|) . Similarly, the outer boundary terms in E^ + E^ (i.e. 
the sums over i) can be controlled on each domain separately. We need 



u\Ay - 2T^)u^ < 0. 
We can therefore take, as in the one-dimensional case, 

P'' = P' = (A+ + 5t)P: , 

where now P^ are projectors to the spaces of eigenvectors of A^ of eigenvalues A+. With these choices the final 
expression for the time derivative of the energy is 



3<0 



(A+ + 2e)ii<-' - «rip - 250||u[, - 



(25) 
(26) 



i>0 



j>0 



and, again, the possible ranges for the different (5's are as in the Id case, Eas. ()9ll0lll|l . 

Notice that nothing special has to be done at a corner, as each direction is treated and controlled independently. 
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C. The general case 

The general case follows the same rules. Namely, we must add penalty terms on the characteristic modes corre- 
sponding to each of the boundary matrices separately and accordingly. 

For example, what to do at the vertices of three patches meeting in the cubed-sphere case discussed later in this 
paper? As we will see, there we have three meshes with coordinates (at a constant radius) {a-^,b^), (a^,6^), and 
{a^,b^) arranged in a clockwise distribution according to the indices, and intersecting at a point. In that case, the 
contribution to the energy (without the penalty terms added to the evolution equations) is proportional to 

{ul^, (A^^ + A''')uIj,) + (4,, + A'^')ul,) + {uIj,, {A^' + A'''')ul^) 

Since the interfaces are aligned to the grids we know that the normals coincide on both sides, therefore we have: 

A"' = A""' A''' = A''' A''' = A""' 

So we include penalty terms on each side, including the end-points of the grids in each direction, (which constitute 
vertices and edges). Note that the characteristic modes at these points are computed with the normal with respect 
to the side that contains this direction. Consequently, points at edges/vertices of a (topologically) cubical grid will 
have two/three penalty terms. 

III. HIGH ORDER DIFFERENCE OPERATORS WITH DIAGONAL NORMS 

In this section we analyze some aspects of Strand's Id difference operators satisfying SBP with respect to diagonal 
metrics, when used in conjuction with the penalty technique to construct high order schemes for handling domains 
with interfaces. 

In particular, we discuss operators with accuracy of order two, four, six and eight at interior points. The requirement 
that these operators satisfy the SBP property with respect to diagonal norms implies that their respective accuracy 
order at and close to boundaries is one, two, three and four, respectively. We will therefore refer to these operators as 
D2-1, D4-2, and i^8-4- Some of these operators are not unique, as the accuracy order and SBP requirements 

still leave in some cases additional freedom in their construction. Indeed, while the first two operators (£'2-1, £'4-2) 
are unique, the De-s one comprises a mono-parametric family, and I?8-4 a three-parametric one. This freedom can 
be exploited for several purposes. For instance, to minimize the operator's bandwidth or its spectral radius. While 
the former produces operators which are more compact, the latter can have a significant impact on the CFL limit 
when dealing with evolution equations. Indeed, for the D^-^ case, minimizing its bandwidth leads to a considerably 
larger spectral radius (though this does not happen in the I?6-3 case) which, in turns, requires one to employ a rather 
small CFL factor for the fully discrete scheme to be stable. 

To analyze this in each case, we numerically solve and discuss the eigenvalues of the amplification matrix of the 
advection equation with speed one, ut — u' , under periodic boundary conditions. The periodicity is imposed through 
an interface with penalty terms and hence the scheme does depend on the penalty parameter S and so will the discrete 
eigenvalues obtained. As discussed in Sectional in the case in which S = —1/2, SBP holds across the interface, 
and the energy for this model is strictly conserved. In other words, the amplification matrix is anti-symmetric and 
the eigenvalues are purely imaginary (see Section^. On the other hand, if (5 > —1/2 there is a negative definite 
interface term left after SBP, and a negative real component in the eigenvalues must appear in the spectrum of the 
amplification matrix (see Section ITTll. 

We additionally discuss the global convergence factor for these operators, and recall a feature associated with the 
mode with highest possible group speed at a given number of gridpoints. Namely, that it travels in the "wrong" 
direction, and that the absolute value of its speed increases considerably with the order of the operator. 

Appendix 1X1 lists, for completeness, some typos in Ref. in some of the coefficients for these high order operators. 

A. Spectrum 

In the following we discuss the range of discrete eigenvalues obtained for the different derivative operators and their 
dependence on S. We pay particular attention on the impact different values of d and the chosen derivative operator 
have on the CFL limit. 
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1. Second order in the interior, first order at boundaries (D2~i scheme) 

Figure 121 shows the eigenvalues obtained using 20, 60, 100 gridpoints and penalty term 5 = — 1/2 (that is, the purely 
imaginary spectrum case, see Section^) for the D2-1 case. The maximum and minimum values are, approximately, 
±1.414, and they seem to be related to the operator near the boundary, as their absolute value does not seem to 
increase with the number of points (instead, the region between the maximum and minimum is filled out). As discussed 
below, in the higher order cases the maximum eigenvalues also seem to be related to the operator near the boundary. 



0.5- 







-0.5 



0.5- 



0- 



-0.5 



0.5- 



0- 



-0.5- 



-1e-10 -5e-11 



5e-11 1e-10 



-1e-10 -5e-11 



5e-11 1e-10 



-1e-10 -5e-11 



5e-11 1e-10 



FIG. 2: Numerically obtained eigenvalues (in the complex plane) corresponding to the D2-1 operator for 5 = —1/2 (purely 
imaginary case). From left to right the plots illustrate the results obtained with a grid containing 20, 60, 100 points respectively. 
It is clear from the figures that these correspond, indeed, to a purely imaginary case. 



Figure 13 in turn, shows the eigenvalues computed with 100 points and 5 = 0,1/10,1/2. A negative real part 
appears, as it should (based on the energy calculation), and the maximum in the imaginary axis slightly decreases 
(to approximately 0.999, not varying much among these three values). However, the maximum absolute value in the 
negative real axis grows quite fast with 5. For example, for 5 = 1/10 such maximum already dominates over the 
maximum in the imaginary axis. The higher order operators analyzed below behave similarly. 



2. Fourth order in the interior, second order at and close to boundaries (7?4_2 scheme) 

Figure 0] shows the equivalent of Figure |3 but now for the £'4_2 case. The maximum is slightly larger than the 
corresponding one for the D2-1 case: approximately 1.936. 

Figure|Sl in turn, shows the equivalent of Figure|31for the current case. As before, a negative real part appears and 
the maximum in the imaginary axis slightly decreases (in this case to roughly 1.371, not changing much among these 
three values of <5). 



3. Sixth order in the interior, third order at and close to boundaries, minimum bandwidth case (Des scheme) 

Figureiniillustrates the equivalent of Figures H2I4|I for the i^g-a case. The maximum is roughly 2.129, slightly larger 
than those of the previous two cases. The behavior for larger values of S is similar to that one found in the previous 
two cases, as seen in FigureEl The maximum in the imaginary axis again decreases slightly compared to the i5 = —1/2 
case (to roughly 1.585) and does not change much among these three values of S. 
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FIG. 3: Eigenvalues corresponding to the D2-1 operator, obtained with a grid containing 100 points. From left to right the 
plots illustrate the behavior for 5 = 0, 1/10, 1/2 respectively. As S becomes larger, a larger (in magnitude) negative eigenvalue 
on the real axis is observed (notice the left-most diamond at y ~ —1.6, —3.5 on the middle and right plots, respectively). 



0- 



0- 



-1- 



—2 J , ^ , , , , ^ , , , , Y . . . . I . . . . I —2 - I ,,,,,,,,,,, Y ■■■■ I ■■■■ I —2 - I ,,,,,,,,,,, Y ■■■■ I ■ 
-1e-10 -5e-11 5e-11 1e-10 -1e-10 -5e-11 5e-11 1e-10 -1e-10 -5e-11 5e-11 1e-10 



FIG. 4: Numerically obtained eigenvalues corresponding to the D4-2 operator for 5 = —1/2 (purely imaginary case). From 
left to right the plots illustrate the results obtained with a grid containing 20, 60, 100 points respectively. It is clear from the 
figures that these correspond, as in Fig|5| to a purely imaginary case. 



4- Eight order in the interior, fourth order at and close to boundaries (Ds~4 scheme) 



The operator has three free parameters, denoted as xi,X2,X3 both in Ref.|y| and here. As mentioned, these 

parameters can be freely chosen to satisfy a given criteria. For instance, they can be fixed so as to minimize the width 
of the derivative operator or yield as small a spectral radius as possible. As we discuss next, these options can yield 
operators with significantly different stability requirements as dictated by the CFL condition. 

Minimum bandwidth operator. The minimum bandwidth case corresponds to the choice (see 0) 



xi = 



1714837 
4354560' 



X2 



1022551 
30481920' 



X3 = 



6445687 
8709120' 



(27) 



Figure El shows for this minimum bandwidth case the eigenvalues for 20, 60, 100 points, for d — 



— 1/2 (purely 
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FIG. 5: Eigenvalues corresponding to the D4-2 operator, obtained with a grid containing 100 points. From left to right the 
plots illustrate the behavior for 5 = 0, 1/10, 1/2 respectively. As S becomes larger, a larger (in magnitude) negative eigenvalue 
on the real axis is observed (notice the left-most diamond at y ~ —2.5, —5 on the middle and right plots, respectively). 
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FIG. 6: Numerically obtained eigenvalues corresponding to the minimum bandwidth De-a operator for S = —1/2 (purely 
imaginary case). From left to right the plots illustrate the results obtained with a grid containing 20, 60, 100 points respectively. 
As in Figs. 12141 . these correspond to a purely imaginary case. 



imaginary case). While for the previous operators we have seen that the maximum eigenvalue increases slightly with 
the order of the operator, in this case the increase is quite large: the maximum is roughly 16.04. This translates into 
a CFL limit for this operator being almost an order of magnitude smaller than the limits for the previous operators. 
Additionally, variation of S does not significantly affect this behavior, as shown in Figure El That figure shows the 
eigenvalues computed with 100 points, and 6 = 0, 1/10, 1/2. The qualitative behavior when increasing d is similar to 
that one of the previous cases. A negative real part appears in the spectrum, and the maximum in the imaginary 
axis slightly decreases, to roughly 16.02, not varying much among these three illustrative values of S. Such a large 
spectral radius for this operator motivates the search for another one, with a more convenient radius at the expense 
of not having the minimum possible bandwidth. 

Optimized operator. We here construct an "optimized D8-4 operator" (which we shall use from here on in the 
D8-4 case) in the sense that it has a spectral radius considerably smaller than that one defined by Ea. H27|) . More 
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FIG. 7: Eigenvalues corresponding to the De-a operator obtained with a grid containing 100 points. From left to right the 
plots illustrate the behavior for 5 = 0, 1/10, 1/2 respectively. As S becomes larger, a larger (in magnitude) negative eigenvalue 
on the real axis is observed (notice the left-most diamond at y ~ —2.5, —6 on the middle and right plots, respectively). 
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FIG. 8: Numerically obtained eigenvalues corresponding to the minimum bandwidth Ds-i operator for 5 — —1/2 (purely 
imaginary case). From left to right the plots illustrate the results obtained with a grid containing 20, 60, 100 points, respectively. 
Although purely imaginary, the maximum (absolute) value in the vertical axis is approximately 16. 



precisely, through a numerical search in the three-parameter space we have found that the following values 

a;i= 0.541, = -0.0675, = 0.748 , (28) 
yield an operator whose maximum absolute eigenvalue in the purely imaginary case (6 = —1/2) is 

This maximum eigenvalue appears to be quite sensitive on these parameters. For example, truncating the above 
values to two significant digits, 

xi = 0.54, X2 = -0.067, X3 = 0.75 , 
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FIG. 9: Eigenvalues for the minimum bandwidth D8-4 operator with S — 0, 1/10, 1/2 (from left to right) and 100 points. Values 
of 5 larger than —1/2 introduce a negative real part in the spectrum but they have little effect on the maximum absolute value, 
which remains at approximately 16. 



gives Xmax = 2.698 and truncating even more, to just one digit, 

xi = 0.5, X2 = -0.07, X3 = 0.7 , 

gives the large value Xmax = 71.76. On the other hand, refining in the parameter search the values in Ea. (|28|l in one 
more digit did not change the maximum of Ea. (|29|l in its four digits here shown. The eigenvalues for 6 = —1/2 for 
this optimized -D8-4 operator, given by the parameters of Ea. (|28|l . are shown in Figure ITUl while Figure ITTI shows 
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FIG. 10: Eigenvalues for the optimized Ds-4 operator with 5 = — 1/2 (purely imaginary case), and 20, 60, 100 points (from left 
to right). Clearly, this modified operator has a much smaller spectral radius, compared to the minimum bandwidth one. 



them for 6 = 0, 1/10, 1/2 and 100 points. As before, a negative real component appears and the maximum in the 
imaginary axis decreases (to around 1.754). 

While completing this work we became aware of similar work by Svard, Mattson and Nordstrom 14], who construct 
an optimized operator with different parameters by minimizing the spectral radius of the derivative itself (rather than 
that of the amplification matrix of a toy problem with an interface, as in our case), obtaining xi = 0.649, a;2 = 
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FIG. 11: Eigenvalues for the optimized Ds-i operator, with 100 points and 5 — 0, 1/10, 1/2 (from left to right). 



—0.104, a;3 = 0.755. When using these parameters in our toy problem with an interface, the resulting spectral radius 
(for twenty gridpoints) is Xmax = 2.241259, while for the parameters we chose [cf. Ea. (|28l) ] is X^ax — 2.241612 )15| . 

B. Global convergence rate 

In general, the global (say, in an L2 norm) convergence factor for these operators will be dominated by the lower 
order at and close to boundaries. However, it is sometimes found that roundoff values for the error in such a global 
norm are reached before this happens, and the convergence factor is different from the one expected from the boundary 
terms. The precise value is found to actually depend on the function being differentiated and whether one reaches 
round-off level. To illustrate the expected behavior in a generic case, we consider the function sin(10a;) + cos(10a;) in 
the domain x £ [0, 27r]. Figure [T^ shows the error (with respect to the exact solution) when computing the discrete 
derivative versus the number of gridpoints, for the difference operators D2-1, D4-2, Dq-3, and Ds-4. The errors in 
the L2 norm are shown until roundoff values are reached (further increasing the number of gridpoints causes the error 
to grow with the number of points involved). Figure [T51 in turn, shows the obtained convergence factors. 

C. Group Speed 

We now turn our attention to the group speed that different discrete modes have when the above considered 
operators are used. To simplify the discussion, we actually restrict ourselves to the periodic case, which lends itself 
for a clean analytical calculation. In this case the operators of order two, four, six and eight satisfying SEP are the 
standard, centered ones {Do, denote the standard centered second order, and forward and backward first 



order operators, respectively): 

15(2) = Do (30) 

= Do{I -h^/QD+D^) (31) 

D^^^ = Doil -h^/6D+D^+h^/30DlDl) (32) 

D^^^ = Dq{I ~ h^/QD+D^ + h^/iQD\D'i - /IAQD\D^_) (33) 

(34) 

In discrete Fourier space, the eigenvalues for these operators are, respectively, 

\2 - sin(C)/C, (35) 

A4 = sin(C)/C(l + 2sin(C/2)/3), (36) 
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Errors with respect to u=sin(10 x) + cos(10 x) 




FIG. 12: L2 norms of the errors obtained when taking the discrete derivative of sin(10a;) +cos(10a;) and comparing it with the 
analytical answer, using 1)2-1, -D4-2, -Dg-s, Dg_4 operators, versus the number of gridpoints. 



Convergence factor with respect to u=sin(1 x) + cos(1 x) 




log^g(Npts) 

FIG. 13: Convergence factors for the curves of Figure [T2l As in this case, the convergence in general will be dictated by the 
lower order the derivative operators have at and close to the boundaries. Note that the lines corresponding to the different 
operators terminate at sequentially fewer points. This is due to the corresponding errors reaching round-off levels, after which 
the convergence factor calculation ceases to have a sensible meaning. 



Ae = sin(C)/C(l + 2sin(C/2)/3 + 8sin4(C/2)/15), (37) 
As = sin(C)/C(l + 2sin(C/2)/3 + 8sin'^(C/2)/15 + 64sin^(C/2)/140). (38) 

where C = Luh, and uj is the associated wave number. 

The highest possible frequency is w = N/2 (with N the number of gridpoints). For that frequency, C — tt and the 
above eigenvalues are all zero. Therefore the mode with highest possible frequency for a given number of points does 
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not propagate. Furthermore, if one examines the group speedy 

d(Xuj) 



(39) 



one finds that at ( = tt this speed is 



= -1, (40) 

= -5/3 « -1.6, (41) 

1,6 ^ -33/15 w -2.2, (42) 

-341/135 w -2.6. (43) 
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Thus, for higher (than two) order operators the velocity of this mode is higher than the continuum one (which is 1). 
But, more importantly, in all cases the speed has the opposite sign. Of course, this effect goes away with resolution, 
since the highest possible frequency moves to larger values as resolution is increased. But, still, is an effect to be taken 
into account. For instance, if noise is produced at an interface, it propagates backwards, and with higher speed. Even 
though this effect is typically very small, it might be noticeable in highly accurate simulations, or in simulations in 
which the solution itself decays to very small values (see Section IIV C 2|l . This could also be a source of difficulties 
in the presence of black holes -or for this matter any system where some propagation speed changes sign- since the 
event horizon traps these high frequency modes in a very narrow region and then releases them as low frequency ones. 
We have observed this in some highly resolved one dimensional simulations, and explains an observed convergence 
drop which goes away when numerical dissipation is turned on. 



IV. TESTS 



In this section we illustrate the behavior of the aforementioned penalty technique, together with the choice of 
different derivative operators. We present tests in one, two and three dimensions. In particular, we implement 
the linearized Einstein equations (off a 'gauge- wave' spacetime |0, 0|) and propagation of scalar fields in black 
hole backgrounds. The former is cast in a way which yields a one-dimensional symmetric hyperbolic system with 
coefficients depending both in space and time while the latter provides an hyperbolic system of equations with space 
varying coefficients and sets a conforming grid for spherical black hole excision. 

Throughout this paper we employ a fourth order accurate Runge-Kutta time integrator. In a number of tests aimed 
at examining the behavior of high order operators we adopt a sufficiently small time step /St so that the time integrator 
does not play a role. Thus, we either choose a suitably small CFL factor or we scale the time step quadratically with 
the gridspacing h. 



A. One dimensional simulations: linearizations around a gauge wave 

As a first test we evolve Einstein's equations in one dimension, linearized around a background given by 

ds^ = e'^^'"('^(^-*»(-dt2 + dx^) + dy^ + dz^. (44) 

This background describes flat spacetime with a sinusoidal coordinate dependence, of amplitude A, along the x 
direction. One of the interesting features of this testbed is that while a linear problem, the coefficients in the 
equations to solve are not only space but also time dependent. 
The non-trivial variables for this metric are 

K.. = :|^cos(7r(a;-t))e^/2^'"(^("-*» , (46) 
a = e^/^^i'^M^-t)) ^ (47) 
= 0. (48) 

We evolve the linearized Einstein equations using the symmetric hyperbolic formulation presented in Ref. [l^ with a 
dynamical lapse given by the homogeneous time-harmonic condition (defined by requiring Di = 0). The formulation is 
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cast in first order form by introducing the variables Ax ■= dxCa/a and dxxx '■= dxQxx- The equations determining the 
dynamics of the (first order) perturbations, which we assume to depend solely on {t,x), are obtained by considering 
linear deviations of a background metric given by Eq. 144|l . That is, we consider 



:/xx , 



9xx 9xx ^ 

^^xx ^^xx I '^^^xx^ 

^xxx ^xxx ^^xxx: 

a = a + Sa, 

A-x — A-x 4" ^Ax , 

replace these expressions in Einstein's equations and retain only first order terms. The resulting equations are 
(henceforth dropping the 6 notation) 

a = -ATTCOs{(t>)a - Kxx + 77^ cos{(j))gxx , (49) 
1 Att 

Ax = - — dxKxx + TTT COs{(j))Kxx , 

a Za 
AtP- 

~ 2d2" + sin(0)) gxx 

+ ^T^COS{(p)dxxx 

Att jAtt'^ 

— cos((/)) A + -Trr- sin(0)a , (50) 

Z Za 

gxx = -AaTrcos{(l))a-2aKxx, (51) 
Kxx = -adxAx - cos{(t>)Ax 
Att"^ 

— (-2 sin(0) + A cos((/i)^) a 

-Att cos{(j))Kxx 
Att 

+ —cofi{(t>)dxxx , (52) 
4a 

dxxx = AaiP' (sin(0) — Acos(0)^) a 

— Aa-K co^{'i))Kxx — Aa'^n coa{(j))Ax 

-2adxKxx , (53) 

where we have defined (p := t:{x — t). This system is symmetric hyperbolic and the symmetrizer used to define the 
energy can be chosen so that the characteristic speeds which play a role in the energy estimate are and ±1. 

We consider here a periodic initial boundary value problem on the domain x £ [—1/2, 3/2], where periodic boundary 
conditions at a: = ^V^i 3/2 are implemented through an interface with penalty terms, as described in Section ITll 

The system must satisfy two non-trivial constraint equations, corresponding to the definition of the variables dxxx 
and Ax (the linearized physical constraints are automatically satisfied by the considered ansatz). When linearized, 
these constraints are 

Q = Cx = -dxgxx + dxxx , (54) 

1 / An \ 
= C'a = Ax - -T {dxa ^cos(0)aj. (55) 

In the first series of simulations we adopt a CFL factor A — 10^'^_39j and consider relatively short evolutions 
corresponding to four crossing times. The -D8-4 derivative is used, and dissipation is added through the dissipative 
operator constructed from —ah7 D'^D'^, suitably modified at boundaries as explained in Appendix IbI so as to make 
it non-positive definite with respect to the appropriate scalar product. Thus, the use of this dissipative operator does 
decrease the order of the spatial discretization by one. The dissipation parameter used is cr = 5 x 10^^. Figure T\M 
exemplifies the behavior observed in the convergence of the field Kxx (the other fields behave similarly). As time 
progresses the convergence order obtained oscillates in a way that is consistent with the accuracy obtained at interior 
and boundary points. 
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time 



FIG. 14: Evolutions of Id linearized Einstein's equations in a periodic domain, with periodicity enforced through an interface 
with penalty terms. Shown is the convergence factor for K^x when using the -Dg_4 derivative, CFL factor A = 10~'^, and 
dissipation cr = 5 x 10"''. While the convergence factors obtained with 41 to 321 points oscillate between the expected order 
at the boundary and that one at the interior, the ones calculated with 321 to 1281 points are not meaningful when the pulse 
is located at interior points as round-off level is reached. 



Next, we adopt as a starting value for the CFL factor defined at the coarsest grip to be A = 0.2 but in subsequent 
grids (refined by a factor of 2) we adopt A = 0.2/2" with (n = 1..3). Figure illustrates the behavior observed; 
again, as time progresses the convergence order obtained oscillates in between the order of accuracy of interior and 
boundary points, with the additional effect of accuracy loss due to the accumulation of error as time progresses. 
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time 

FIG. 15: Same as previous figure, but with a decreasing CFL factor given by A = 0.2/2" with (n = 1..3). 



Finally, we compare the above results with those obtained in the "truly periodic case" , ie. when periodicity is used 
explicitly to employ the same derivative operator at all points. We again consider cases where a sufficiently small 
CFL factor (= 10"'^) is used or the time-step is scaled quadratically. Figures [TBI illustrates the observed convergence 
rate for the field Qxx- As above, dissipation is added through a seventh order dissipative operator (but now with 
no modification at boundaries needed) with same dissipative parameter: cr = 5 x 10""*. While the errors remain 
above round-off level the observed convergence rate is consistent with the expected one of seven, as the orders of the 
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derivative and dissipative operators employed are eight and seven respectively. Certainly, dissipation of higher order 
could have been introduced by simply employing the KO style operator h^{D_D+)^ , but we have adopted this one to 
more directly compare with the case with interface boundaries. For the highest resolutions the errors reach round-off 
level and the obtained convergence factors yield non-sensible values. 

Convergence factor (D ) Convergence factors (D ) 



Gauge wave, periodic boundary conditions, >i=0.00l Gauge wave, periodic boundary conditions, quadratic X 




time time 



FIG. 16: This figure shows evolutions similar to those of Figs ll4l 1151 with the only difference that periodicity is here enforced 
explicitly. The convergence factors for the metric component g^x are shown. 



As an illustration of what is observed with other derivatives, we briefly discuss some simulations using the I?6-3 
operator and a fixed CFL factor (given, as before, by A = 10"'^). Analogously as to was done above, dissipation is 
here added by extending -as discussed in Appendix IbI- the operator ah^D^D\ at and near boundaries in order to 
make it non-positive definite with respect to the appropriate scalar product; a dissipation parameter a = 10~^ is 
used. The observed results are illustrated in Figure [T7I which shows the self-convergence factor for K^x- As before, 
it oscillates between the order of the scheme in the interior and that one at boundary points. 




FIG. 17: Similar to Fig. 1141 but using the De-s derivative. The convergence factor for K^x is shown. 



Summarizing, the results presented indicate that, in the case where boundaries are present, the worst case scenario 
-as far as the expected convergence rate relates- is determined by the accuracy order at and close to boundary points. 
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Two and three-dimensional simulations. Problem set-up 



In this section we solve the wave equation for a scalar field (f) propagating on a fixed background, 

V^Va^^O, 

where V is the covariant derivative associated with the metric of the background. We will consider two backgrounds: 
a 2d one consisting of the unit sphere with its standard metric and a 3d one consisting of a rotating Kerr black hole 
background. 

We start by describing in more detail the equations solved and the multiple coordinate system used, and then 
present the actual results of the simulations. 

A strictly stable scheme for the wave equation in a time-independent, curved spacetime 

The wave equation in a time-independent background can be written, using any coordinates in which the metric is 
manifestly time independent as, 

^ ^ an (56) 
n = l3'a-^D,{an) + h-^/'^D,{h^/^l3m + ah^/^H'^dj) (57) 
d, = A(a7r) (58) 

where H^^ := h^^ — a~^/3*/3^ , h^^ is the inverse of the three-metric, h = det(/i.y ), a is the lapse and /3* the shift vector. 
The advantage of writing the equations in this way is that one can show that if D is any difference operator satisfying 
SBP, this form of the equations guarantees that the semidiscrete version of the physical energy is a non-increasing 
function of time. When the killing field is timelike this means that there is a norm in which the solution is bounded for 
all times, thus suppressing artificial fast growing- modes without the need of artificial dissipation (see 1^ for details). 

We now look at the characteristic variables and characteristic speeds with respect to a "coordinate" observer. That 
is, the eigenfields and eigenvalues of the symbol A^rii, where denotes the principal part of the evolution equations 
and the normal to the boundary liSI. The characteristic variables with non-zero speeds 

= (±a + f3''nk)ih'^n,nj)^/^ 

[where fik — nk{h'^-' ninj)^^^^] are 

v"^ = X'^U + aW^h.dj ; 

while the zero speed modes are 

v'^ = di — fiidjiiP . 



Cubed-sphere coordinates. 

The topology of the computational domain in our 2d simulations is 5^, the unit sphere, while in our 3d ones it is 
S'^ X . Since it is not possible to cover the sphere with a single system of coordinates which is regular everywhere, 
we employ multiple patches to cover it. A convenient set of patches is defined by the cubed sphere coordinates^ defined 
as follows (for a related definition see for instance [20l|). 

Each patch uses coordinates a, 6, c, where c = a/x^ + + z^, the standard radial coordinate, is the same for the 
six patches (x,y,z are standard Cartesian coordinates). The other two coordinates, a, 6 are defined as 

• Patch (neighborhood ofa; = l): a — z/ x,b — y / x 

• Patch 1 (neighborhood of y = 1): a ~ z/y^ b — —x/y 

• Patch 2 (neighborhood of x = —1): a = —z/x, b = y/x 

• Patch 3 (neighborhood of y = —1): a = ^z/y, b — —x/y 

• Patch 4 (neighborhood of z = 1): a — —x/z, b = y/z 
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• Patch 5 (neighborhood of z = —1): a = —xjz, b = —y/z 
Similarly, the inverse transformation is: 

• Patch 0: x = c/ D,y = cb/D, z = ac/D. 

• Patch 1: X = —bc/D, y = c/ D, z = ac/D 

• Patch 2: X — —c/D, y = —cb/D, z = ac/D. 

• Patch 3: x — be/ D,y — —c/ D, z = ac/ D 

• Patch 4: X — —ac/D, y — cb/ D, z — c/ D 

• Patch 5: X — ac/D, y — cb/D, z — —c/D 

with D :— \/l + a^+l?. This provides a relatively simple multi-block structure for which can be exploited to 
implement the penalty technique in a straightforward manner. Each patch is discretized with a uniform grid in the 
coordinates a and b, and the requirement of boundary points coinciding in neighboring grids is indeed satisfied. Figure 
IT51 shows this gridstructure, for 20 x 20 points on each patch. 

1 



0.5 



z 



-0.5 



-1 



FIG. 18: Cubed-sphere coordinates for . 



B. 2d Simulations 

We now discuss simulations of the wave equation on the unit sphere in cubed-sphere coordinates, written in strictly 
stable first order form [Eas. 1)56157158(1 ]. The metric used, therefore, is flat spacetime projected to the r = 1 slice, 
which in local coordinates is 

ds^ = -dt^ + D^'^ [(1 + b^)da'^ + {l + a^)db'^ -2ab da db] , 

where D := ^/(l + + 62), 

Figure El shows simulations using the D4-2 derivative and its associated dissipative operator constructed in Ap- 
pendix^ which we call K06, using n x n points on each of the six patches, where n — 41,81, 161,321. The initial 
data for IT corresponds to a pure / = 2, m = 1 spherical harmonic. The CFL factor used is A = 0.125 and for each set 
of runs two values of dissipation are used: a = 10~^ and a = 10~^. As can be seen from the Figure, the self conver- 
gence factor obtained with these resolutions is above the lower value (two) expected from the order at the interfaces. 
The reason for the lower order at the interfaces not dominating is likely due to the fact that the initial data is an 
eigenmode of the Laplacian operator, and the solution at the continuum is just an oscillation in time of this initial 
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data, without propagation across the interfaces. Indeed, the oscillations in the convergence factors in Figll9l appear 
when the numerical solution goes through zero, and the frequency at which this happens coincides approximately 
with the expected frequency at the continuum for this mode. 

The same initial data is now evolved with the De-s derivative and K08 dissipation (again, see Appendix and 
the results are shown in Figure EOl As before, A = 0.125 and n = 41,81, 161,321 points are used, but the values of 
dissipation shown are now a — (i.e., no dissipation) and a = 10^'^. At the same resolutions there is a small difference 
in the obtained convergence factors, depending on the value of ti, but with both values of this parameter the order 
of convergence is higher than the lower one expected from the time integrator if this one dominated. Figure I2UI also 
presents a comparison made with a smaller CFL factor: A — 0.0125, keeping the dissipation at a — 10^^. One more 
resolution is used (641 points) to look for differences between the solutions obtained with the two CFL factors, but 
they do not appear. This seems to suggest that at least in this case, and for these resolutions, it is not necessary to 
use too small a CFL factor in order to avoid the time integrator's lowest order to dominate over the higher spatial 
discretization (see Fig l23l for another instance where this happens). It is also worth pointing out that the difference 
between the two highest resolutions is not quite at roundoff level, but it is rather small (of the order of 10^^ if scaled 
by the amplitude of the initial data), as shown in Figure 1^ That figure shows the L2 norm of the differences between 
the solution at different resolutions, for the simulations of Fig. I^Olwith A = 0.0125 and a = 10~^. 

2D Convergence Factor D4-2/K06 
1=2, m=1 initial data 

^1 ' ^ \ 5 r 1 

41-81-161, a=10 , X=0.125 

— 81-161-321, 0=1 0"^ A, =0.1 25 
"■^^ — 41-81-161, a=10'^, X=0.125 

— 81-161-321, 0=10"^ ?i=0.125 




2.5 - 



time 

FIG. 19: Evolutions of the wave equation on the unit sphere, using cubed-sphere coordinates and a pure I — 2,m = 1 spherical 
harmonic as initial data. Shown is the convergence factor when the D4-2 derivative and the K06 dissipation operators are 
used. 



Finally, Figure [23 shows evolutions of the same initial data, with the Ds~4 derivative, and no dissipation. The CFL 
factor is decreased when resolution is increased, much as in Section llV Al so that the order of the time integration 
does not dominate over the higher one of the spatial discretization. That is, for the resolutions shown we used 
A = 0.25, 0.125, 0.0625, 0.03125, 0.015625. The convergence factor obtained is also higher than that one expected from 
the lower order at the interfaces (presumably for the same reason as before, the initial data adopted) and higher than 
that one of the time integration. 

C. Three dimensional simulations 

In this application we consider fields propagating on a Kerr black hole background, as governed by equations 
(|56I57I58|I . with the background metric written in Kerr-Schild form and cubed-sphere coordinates used for the an- 
gular directions. Homogeneous maximally dissipative boundary conditions are used at the outer boundary, while 
no condition is needed at the inner one if it is appropriately placed inside the black hole so that it constitutes a 
purely-outflow surface. 
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2D Convergence Factor D6-3/K08 
1=2, m=1 initial data, cfl=0. 125 



2D Convergence Factor D6-3/K08 
l=2,nn=1 initial data, a=10 ^ 
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FIG. 20: Same evolution equation and initial data as those used in Figure IT??1 except that now the De-s derivative and K08 
dissipation are used. 



Difference in n at different resolutions, D6-3/K08 
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FIG. 21: 1/2 norms of the differences between the solution at different resolutions, for the simulations of Figure 1201 with 
A = 0.0125 and a = 10"^ 



The Kerr-Schild metric in cubed-sphere coordinates 
The Kerr metric in Kerr-Schild form is 

ds'^ = rjf,^ + 2Hl^lydx^dx'' 
where -q^i, is the flat metric [with signature (— , +, +, +)], 



H = 



+ cos 



„2 l/„2 a2\ '1 



(59) 
(60) 
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2D Convergence Factor D8-4 (no dissipation) 
1=2, m=1 initial data 
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FIG. 22: Simulations of the same equation and initial data as those of the previous figures, except that now the higher order, 
D8-4 derivative is used and no dissipation has been added. 



= x^+y' + z\ (61) 

and Z'^ is a null vector (both with respect to the flat metric and whole metric). 

Therefore, in order to write the above metric in cubed-sphere coordinates one needs to write ry^i^ and Z^ in these 
coordinates. A straightforward change of coordinates of the first one gives 

ri^^ = -^dt"^ + dc^ + c^D-'^ [{l + b^)da'^ + {1 + a^)db^ - 2abdadb] 

with D := v/(l + a2+&2). 

In Cartesian coordinates the co-vector is, in turn, 

rx + Ay ry — Ax z 
I = l^dx^' =dt+ dx + dy+-dz 

with x'' — [t, X, y, z) which, when changed to cubed-sphere coordinates [t, a, 6, c) gives 

, -c^aA^ia^ - D'^) , -~c^A(D^r + Aba^)„ c(D^r^ + a'^A^) , , , 
I ^ dtA — — ^ TTTT-^da H j-—: — — -db + ' , ^ T^rr-dc for patches 0-3 (62) 

, -^A{D^rb + aA) ^ c^A{raD^ - Ab) „ c{D^r^ + A^) , , , 

I ^ dt + — . \ „ T^da + ^ ' , „ T^db + , „ -pr^dc patch 4 63) 

-c^ A[-D'^Tb ^ aA) ^ -c'AiraD^ +Ab) c{D^r^ + A^) ^ 

I = dt-\ — da H ^. , „ — db + ^„ , „ T^dc patch 5 (64) 

D'^r{r^ + A^) D^r{r^ + A^) D^r{r^ + A^) ^ ' 

To write the wave equation, one also needs the inverse metric, which is 
where all indices are raised with rf"^ (the inverse of the flat metric). The non-zero components of the latter are: 

.r^^^^, (65) 

.^^^^M, (66) 

7?- = 1, (67) 
77" - -1, (68) 

= (69) 
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and the vector in the cubed-sphere coordinates is, 

~aA{rb-A) ^(a^ - D^) dD^^+a^A^) 



If- 



-1, ■ 



r(r2 + A'^) 



-1, 



-1, 



A{rb + aA) A{ar - bA) c{D^r^ + A^) 

A{rb-aA) A{ar + bA) c{D^r^ + A"^) 
r{r^ + A2) ' " r(r2 + A2) ' D^r{r^ + A^) 



patches to 3 
for patch 4 
for patch 5 



(70) 

(71) 

(72) 
(73) 



1. Convergence tests 

Figure 051 shows the differences, in the L2 norm, between the numerical solutions at consecutive resolutions, using 
the I?8-4 scheme, with no dissipation. The number of points in the angular directions is kept fixed to 16 x 16 points 
on each of the six patches, and the number of radial points ranges from 101 to 6401. The background is defined by a 
non-spinning black hole, and the inner and outer boundaries are at 1.9M and 11. 9M, respectively. Non-trivial initial 
data is given only to IT, in the form of a spherically symmetric Gaussian multipole: 



n(0, x) = Aexp(r - ro)^/CT, 



' 



(74) 



with rp — 5M,a — M,A = 1. The resulting convergence factors, the normalized differences \\un — W2Af||/||'"Jv|| 
and the non- normalized ones, Hu^v — M27v||; are shown. The use of multiple patches not only allows for non-trivial 
geometries, but additionally one is able to define coordinates in a way such that resolution is adapted to the problem 
of interest. For example, in the geometry being considered, x i?, one employs a number of points in the angular 
direction limited by the expected multipoles of interest and concentrates resources to increase the number of points 
in the radial direction. As an example, the relative differences between the solution at different resolutions shown 
in Figure 1221 reaches values close to roundoff, with modest computational resources. Even though the solution here 
evolved is spherically symmetric at the continuum, as discussed below the number of points used on each of the six 
patches that cover the sphere can reasonably resolve an Z = 2 multipole. Next, Figure shows similar plots, but 
keeping the number of radial points fixed (to 101), using Na x Na points on each of the six patches in the sphere, 
with Na = 21,41,81. 



2. Tail runs 



To illustrate the behavior of the described techniques in 3d simulations we examine the propagation of scalar fields on 
a Kerr black hole background. The numerical undertaking of such problem has been previously treated using pseudo- 
spectral methods ^l|, which for smooth solutions allows the construction of very efficient schemes. As explained next, 
the combination of multi-block evolutions with high order schemes also lets one to treat the problem quite efficiently. 
A detailed study of this problem will be presented elsewhere ; we here concentrate on two representative examples 
of what is achievable. 

In the first case we examine the behavior of the scalar field propagating on a background defined by a black hole 
with mass M = 1 and spin parameter a — 0.5. Non trivial initial data is given only to H, with a radial profile given 
by a Gaussian pulse as in Ea. (|74|l and angular dependence given by a pure I = 2 multipole. The inner and outer 
boundaries are placed at r = 1.8M and r — 1001. 8M respectively. We adopt a grid composed by six cubed-sphere 
patches, each of which is discretized with 20 x 20 points in the angular directions and 10001 points in the radial one. 
This translates into a relatively inexpensive calculation. 

We adopt the Ds-4 derivative operator, add no artificial dissipation and choose a CFL factor A = 0.25. The salient 
features of the solution's behavior observed are summarized in figure 1^ which shows the the time derivative of the 
scalar field, as a function of time, at a point in the equatorial plane, on the even horizon. At earlier stages, the familiar 
quasi-normal ringing is observed. Next the late-time behavior of the field reveals the expected tail-behavior as a fit 
in the interval t G [35QM,750M] gives a decay for 11 of which agrees quite well the expected decay of t~*. 

This can be understood in terms of the generation of an / = mode in the solution due to the spin of the black hole 
[2lll2^ . Finally, noise can be observed appearing at i ~ 800Af due to the outer boundary. This noise, however, is not 
related to physical information propagating to the outer boundary and coming back (for this one would have to wait 
till i « 2, OOOM) but, rather, is related to spurious modes with high group velocities traveling in the wrong direction, 
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Difference in n at different resolutions, D8-4 (no dissipation) 
3D simulations, 1=0 initial data 
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Difference in n at different radial resolutions, D8-4 
3D simulations, 1=0 initial data, ?i=0.25, a=0 
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FIG. 23: Convergence test in the radial direction, with CFL factors A = 0.25, 0.025. Notice that no appreciable difference is 
found between A — 0.25 and a smaller value and even with A — 0.25 the convergence factor is not dominated by the time 
integrator. 



as described in Section lTTTl As discussed there, for an eighth order centered derivative the speed of this spurious mode 
is around —2.6, which roughly matches with this noise appearing at i ~ 800Af . We have checked that this boundary 
effect does go away with resolution, by introducing some amount of dissipation or by pushing the outer boundary 
farther out. Notice that while the first two options allow one to observe the tail behavior for much longer, eventually 
physical information would travel back from the outer boundary and "cavity" effects which affect the decay would 
take place. Indeed, the behavior would no longer be determined by a power law tail but by an exponential decay |24| . 

Figure 1201 shows a similar run, in this case however the black is not spinning. A fit to the solution in the tail regime 
gives a decay for (/) of t^^-^^, which again matches quite well the expected decay of psf . 



V. FINAL COMMENTS 



As illustrated in this work, the combination of the penalty technique together with those guaranteeing a stable single 
grid implementation for hyperbolic systems provides a way to achieve stable implementations of multi-block schemes 
of arbitrary high order. A similar penalty technique for multi-block evolutions is also being pursued in conjuction 
with pseudo-spectral methods p^ . 

The flexibility provided by multiple grids can be exploited to address a number of issues currently faced in simula- 
tions of Einstein's equations, among these 

• The desire for a conforming inner boundary. This plays a central role in ensuring a consistent implementation 
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Difference in n at different angular resolutions, D8-4 
3D simulations, 1=0 initial data, ?i=0.25, a=0 
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FIG. 24: Convergence test in the angular direction, with a CFL factor A = 1/4 



Tail Run (a=0.5,l=2) 




time [M] 

FIG. 25: Behavior of the time derivative of the scalar field for a pure multipole Z = 2 initial data, on a Kerr background, 
with a = 0.5. Inner and outer boundaries are at l.SM and l,001.8Af, respectively and the six-patches grid is covered by 
20 X 20 X 10, 001 points on each patch. The Ds-4, derivative is used, with no artificial dissipation. The noise at f ~ 800M is due 
to the mode with group speed —2.6 discussed in Section UTTI hitting the outer boundary and reaching the observer at the black 
hole horizon. This noise goes away by either pushing the outer boundary, increasing resolution and/or adding dissipation. The 
average slope for 11 in the interval t £ [350M, 750A/] gives a decay for 11 of t"*"^, in good agreement with the expected decay 
of The inset shows a zoom in at the tail behavior. 



of the excision technique together with a saving in the computational cost of the implementation. 

• The need for a smooth outer boundary. This removes the presence of corners and edges which have proved 
difhcult to dealt with even at the analytical level ^27, 28] . Furthermore, a smooth outer boundary simplifies 
tremendously the search for an efficient matching strategy to an outside formulation aimed to cover a much 
larger region of the spacetime with a formalism better suited to the asymptotic region (see, for example, p9| 
and |3(ll3lLl3a|). 

• The use of a grid that is better adapted to the description of wave phenomena as they propagate in the region 
far from the sources. 
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Tail run (a=0,l=2) 





400 
time [IVI] 



FIG. 26: A simulation with the same parameters as those of Fig. 1251 but with a non-spinning black hole. The average slope 
for $ in tail regime gives a decay of t~^'^^ , again in good agreement with the expected decay of 
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APPENDIX A: COEFFICIENTS FOR HIGH ORDER OPERATORS WITH DIAGONAL METRICS 

For completeness we point out here some misprints in Ref.^] in some of the expressions for the diagonal metric 
cases. 

• -D2-1: No typos. 

• D4-2'- It says a2 — ^1/2, but it should be 0:2 — —1/12. The expressions for q2i are also missing; they should 
be the following: 

8 59 

920 — — 921 = —7-7 = 923 — —921 , 924 — —920 

86 80 

• Dq^^: No typos in the scalar product or coefficients for the derivative, neither in the general case nor in the 
minimum bandwidth one. 

• Ds-i'. The operator that has the minimum bandwidth is correct, but the three-parametric one has a typo in 
one of the coefficients (the scalar product is correct): it says 

906 = 49(-1244160a;i + 18661400x3 - 13322233)/17977668 
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when it should say 

qo6 = 49(-1244160xi + 18662400x3 - 13322233)/17977668 

APPENDIX B: DISSIPATION FOR HIGH ORDER DIFFERENCE OPERATORS WITH DIAGONAL 

NORMS 

The addition of artificial dissipation typically involves considering a dissipative operator Qd which is non-positive 
with respect to the scalar product with respect to which SBP holds, i.e., 

<u,QdU><Q Vu. (Bl) 

We start the construction of such operators satisfying this property by defining them as 

Qd = i-ir-'ah'^iD+D-r (B2) 

on gridpoints lying within the range [r, s] contained in the interval in which the weight used in the scalar product is 
one. That is, the interval in which the difference operator is one of the centered ones of Eas. (|30l31l32l33|) : for example, 
if the gridpoints range from to N, then (r, s) must satisfy 

r>l,s<N-l forD2-i; (B3) 
r > 4, s < iV - 4 for Di_2 ] (B4) 
r>6,s<N-6 for Des ] (B5) 

r > 8, s < iV - 8 for Ds-a ■ (B6) 

As we will see later, in some cases our construction of dissipative operators imposes stricter constraints on the range 
of allowed values for r, s for each derivative. 

If n = 2m — 1, the operator (|B2|) is the standard Kreiss-Oliger dissipation (KO) ,33], which is negative definite 
in the absence of boundaries (when the weight in the scalar product is identically one). The choice n = 2m — 1 
ensures that the added dissipation has the same 'scale' as the principal part (that is, length^^) and that the resulting 
amplification factor is independent of resolution. 

For each of the derivatives D2-1, £'4-2, -De-s, -08-4 7 we seek to extend Qd in Ea.HB2|l so that the resulting operator 
is negative definite with respect to the corresponding SBP scalar product. As mentioned above, we denote the points 
in which Qd is given by Ea. (|B2p a,s i = r . . . s. 

The identities of Appendix C are used in the calculations needed for the construction of the operators be- 
low. These identities let one express the norm of the dissipation in the interior, which is proportional to 
(~l)"-i(u, (L>+L)_)"it)['~^"l, as 

r-l N 

(-1)"-1(m, {D+D^)"u)^'''^'^ non-positive definite terms + ^u^(. . .) + ^u^(. . .) . (B7) 

s+l 

Once this is done, the norm of the whole dissipative operator can be written as 

r-l N 

{u, Qrfu)^"'^' = non-positive definite terms -I- Ui[{. . .)i + haiQdUi] + Ui[{. . .)i + haiQdUi] (B8) 

s+l 

and it is straightforward to make it non-positive definite. For example, by choosing 

Qdu^ = , (B9) 

which is how we proceed below. Notice however this is not the only way to proceed. For instance, one could try to 
make the two sums of Ea. (jB8|) cancel (as opposed to requiring the terms cancel at each gridpoint). 

A more general approach for constructing dissipative operators that are negative definite with respect to the 
appropriate (SBP) scalar product has been recently presented by Mattson, Svard and Nordstrom [s^ . 
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(u, Qdu)^'^' = -ah"'\\D+D_u\\f^^j^ ^^^+h'^aiUiQdU.i + h'^(7iU,_Qa 



1. Fourth derivative dissipation, and KO type for D2 i 

In this case the interior operator is 

Qd^-ah"(D+D_)\ 

One can then spht the norm of the dissipative operator over the whole grid into terms involving the left and right 
ranges [0,r], [s,N] (of the yet undefined operator) and the interior terms: 

r-l AT 
s+1 

where cr, are the scalar product weights. Using the identities of Appendix lUl One can see that 

r-3 N 

IdUi 

s+3 

+ hUr-2 {<Jr-2Qd + f7/l""^L'+) U,.-2 + hUr-1 [(Jr-lQd ~ ah"^'^{2D+D- - Z?^]) Ur-1 

+hus+2 {crs+2Qd + cFh'^^'^D'i) Us+2 + hu^+i [<Ts+iQd - <Th''-^{2D+D^ - D^)] u^+i 
There are several options now to control the contribution of these terms, the simplest one is 

QdUr-2 = D']ur-2 

Of -2 

QdUs+i = {2D+D^- Dl)us+i 

QdUs+2 = --^^/l""^L'^Us+2 

and Qd = everywhere else, which implies r > 2, s < iV — 2. Notice that, as anticipated, this is a further constraint 
on the possible range of r, s [c.f. Ea. (|B3|) ]. The operator constructed not only satisfies the non-positivity requirement, 

but also transforms according to D_ — > —1?+ under the symmetry i > N — i (as it should). 

Preferred choice: 

A KO type dissipation for -D2-1, which we call K04 (because it involves fourth derivatives in the interior) corre- 
sponds to 

n = 3 . 

With this choice the order of the scheme is one in at least two points (while in the absence of dissipation would be of 
order one in just one point). In general, the best option for r, s is 

r = 2,s = iV-2, 

as in that case the dissipation operator is non-trivial everywhere (while choosing r>2,s<A^ — 2 would imply that 
the dissipation is zero at -and even possibly close to- the boundaries) . With this choice for r, s the order of the 
scheme is one at the last two boundary points. 

2. Sixth derivative dissipation, and KO type for D4-2 

We now proceed to construct a higher derivative dissipative operator. As before, we adopt at the interior points 
the standard dissipation operator, 

Qd = ah"DlDi 

Then, defining v = D^D^u and employing the identities of Appendix lUl one obtains 
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r-l N 

s+l 
+Ur-3h + (Jr-sQdUr-s) 

+Ur-2h {ah"~^{-hDl_ + 2D+)Vr-2 + (yr-2QdUr-2) 

+Ur-lh [cr/l"~3 {h{2D+D- - Dl)Vr-l - D+Vr-2) + ar-lQdUr^l] + 

+Us+ih [cr/i"-3 [h{2D+D^ - Dl)vs+i + D^Vs+2) + as+iQdUs+i] 

+Us+3h {ah'^^^D^Vs+2 + CTs+sQdWs+s) ■ 
Here again one has several options, the simplest one is: 

QdUr-3 ^ U+Vr-2 

QdUr-2 = [hD\-2D+)Vr-2 

Cr — 2 

QdUr-l = [h{2D+D- - DX)Vr-l - D+Vr-2] 

QdUs+i = [hi2D+D^ - Dl)v,+i + D^Vs+2] 

QdUs+2 = {hD^_+2D-)vs+2 

O's+2 

QdUs+3 = D^Vs+3 

and Qd = everywhere else. Notice that, as in the D2-1 case, Qd transforms under the symmetry i N ~ i Qd as 
it should. Notice, however, that since we are modifying the dissipative operator in three points near the boundary, 
no further constraint occurs in the allowed range for r, s [c.f. Ea. ljB4|l ]. unlike the -D2-1 case. 

Preferred choice: 

A KO type dissipation for D4-2, which we call K06 corresponds to 

n = 5 

In general, the best choice for r, s is 

r = 4, s = 7V-4, 

as in that case the dissipation does not reduce the order of the overall scheme ^Jl- However, a drawback of our 
simplified method for constructing dissipative operators (which cannot be remedied by adopting different values for 
r, s) is that in this case the dissipation is zero at the last gridpoint. 

3. Eighth derivative dissipation, and KO type for De-s 

As usual, we start with 

Qd^ -ah^^D^Di i = r...s (BIO) 



32 



Then, using the identities of appendix O the norm of the dissipative operator results 

r-l N 
s+l 

+Ur-i [-crh"~^ (h^{2D+D_ - D\)wr-i - hD+Wr-2 + -D+a,._2) + crr-ihQdUr-i\ 

+Ur-2 [-ct/i"~^ {{-h^D\ + 2hD+)Wr-2 - 'iD+ar-2) + (^r-2hQdUr~2] 
+Mr-3 [-(jh"'~^{-hD+Wr-2 + 3£'+a.r-2) + ^r-S ^Qd^r-s] 

+Mr-4 {(Th"^^D+ar-2 + (Tr-AhQdUr-i) + Ug+A {crh"-^^ D+as+2 + 0's+4^Q£iWs+4) 

+Ms+3 [-Cr/l""^(/li:'_U;s+2 + 3£'+as+2) + CTs+S^QdMs+s] 

+Us+2 [-Cr/l""^ {~Q?D^_ + 2hD^)Ws+2 - 3D+a^+2) + crs+2/lQdWs+2] 

+its+i [-cr/i""^ {h^{2D+D^ - D'i)ws+i + hD^Ws+2 + £'+as+2) + CT^+i/iQdUs+i] 

where w = {D^D )'^u,a = D^D^D-u. The simplest choice is given by setting to zero all these coefficients, which, 
when expanding w and a gives 

QdUr-l = [h'^{2D+D_ - D\){D+D^fur-i - hD+{D+Djfur-2 + {D+D.fur-2] (BU) 

QdU,_2 = [{-h^Dl + 2hD+){D+D_f - 2,{D+D_f] u,_2 (B12) 

QdUr^3 - [-/ii?+(I?+I?_)4 + 3(Z?+i?_)'] U,_2 (B13) 

QdUr-4 = {D+D^fur-2 (B14) 

O'r-4 

Qd".+4 = (i?+i?-)2u,+2 (B15) 

O's+4 

(j/,"-4 

QdU,+3 = [hD^{D+D^)^ + 3iD+D^)^]us+2 (B16) 

(j/,n-4 

QdM,+2 = [-(/i^i^^ ^ 2/^i^_)(i?+i^_)4 - 3{D+D^)^] Us+2 (B17) 

O-s+2 
(j/jn-4 

aQdW.+i = [/i2(2D+£)_-D^)(D+£>_)4u,+i + ;iD_(D+D_)4u^_^2 + (£'+£'-)^u.+2] (B18) 

Cs+l 

and = everywhere else. 
Preferred choice: 

If we want a KO-type dissipation for the -De-s derivative, which we call K08, we have to choose 

There are six points near each boundary where the difference operator has order three, while the above dissipation 
has four points where that happens. Therefore, the order of the whole scheme is not spoiled if 

r = 6, s = -/V — 6 
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As in the D4-2 case, the drawback of our construction is that the dissipation is zero near boundaries; in this case at 
the last two points. 

4. Dissipation for Ds-i 

In this case we have not been able to write an interior KO dissipation which does not spoil the eighth order accuracy 
of the derivative (that is, one as in Ea. l)B2|l . with m = 5) in the form of Ea. ljB8|l . Therefore our simplified approach 
does not work in this case. However, since the expressions of the previous subsections are valid for a general scalar 
product, we can use, for example, expressions l|B10|l and ljBlllB18|l with the weights corresponding to the -D8-4 
derivative, and 

n = 7,r = 8,s = iV-8 

This results in a KO dissipation of the K08 type which is zero at the last four points, with the drawback that the 
order of the scheme is reduced to seven at the interior points, and three near boundaries (as opposed to eight and 
four, respectively, in the derivative itself). 

APPENDIX C: USEFUL PROPERTIES IN THE CONSTRUCTION OF DISSIPATIVE OPERATORS 

It is straightforward to show that with respect to the scalar product and norm 

the following properties hold: 

{u,D+v)^^''^ = -(I?_u,w)['-+i'^+il (CI) 
iu,D_v)^-''^ = -iD+u,v\,.-i.s-i]+UjV^\:._, (C2) 

{u, Dqv)^'''''^ = -{DoU,v)[r,s]+ ^{'UjVj + l+Uj + iVj)\-'r_i (C3) 

The proofs for the first and third identity can be found in Ref. 35] , and the second one is trivially obtained from the 
first. 

The following equalities are also straightforward to check, though obtaining them becomes increasingly more cum- 
bersome as more derivatives are involved. 

— —{D-U, D-V)^^''^^^^ — Ur-l-D+Wr-l + U^+l-D-Us+l 

(u,(D+L'-)=^z;)['''"l = {D+D-u,D+D-vf-'^-''+^^ + 

+ ^ ({2D+D^ Dl)Wr-l + ^ (-Dl + ^D+) "^D+Wr^, 



{2D+D^ - Di)ws+i + \d^Ws+^ + ^ (-Dl - ] w,+2 + 



(u, (D+i?_)S)['-'^l = iD^q,D-a)[r,s+2] 
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+ ^ (^{2D+D- - Dl)Wr-l - ^D+Wr-2 + ^D+ar-2 

Ur_3 / 3 
H I -D+Wr-2 + ^-D+a,._2 



H ^ I L'-Ws+2 + J^L)+as+2 I 

where w = D+D^v,p = D+D^u^q = D^D+u^a = D^D+v, and denotes the value of the scalar product at 
gridpoint i. 
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